home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / zmlri.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  9.2 KB  |  259 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((zeror 0.0) (zeroi 0.0) (coner 1.0) (conei 0.0))
  12.   (declare (type double-float conei coner zeroi zeror))
  13.   (defun zmlri (zr zi fnu kode n yr yi nz tol)
  14.     (declare (type (simple-array double-float (*)) yr yi)
  15.              (type f2cl-lib:integer4 kode n nz)
  16.              (type double-float zr zi fnu tol))
  17.     (prog ((i 0) (iaz 0) (idum 0) (ifnu 0) (inu 0) (itime 0) (k 0) (kk 0)
  18.            (km 0) (m 0) (ack 0.0) (ak 0.0) (ap 0.0) (at 0.0) (az 0.0) (bk 0.0)
  19.            (cki 0.0) (ckr 0.0) (cnormi 0.0) (cnormr 0.0) (fkap 0.0) (fkk 0.0)
  20.            (flam 0.0) (fnf 0.0) (pti 0.0) (ptr 0.0) (p1i 0.0) (p1r 0.0)
  21.            (p2i 0.0) (p2r 0.0) (raz 0.0) (rho 0.0) (rho2 0.0) (rzi 0.0)
  22.            (rzr 0.0) (scle 0.0) (sti 0.0) (str 0.0) (sumi 0.0) (sumr 0.0)
  23.            (tfnf 0.0) (tst 0.0))
  24.       (declare
  25.        (type double-float tst tfnf sumr sumi str sti scle rzr rzi rho2 rho raz
  26.         p2r p2i p1r p1i ptr pti fnf flam fkk fkap cnormr cnormi ckr cki bk az
  27.         at ap ak ack)
  28.        (type f2cl-lib:integer4 m km kk k itime inu ifnu idum iaz i))
  29.       (setf scle (/ (f2cl-lib:d1mach 1) tol))
  30.       (setf nz 0)
  31.       (setf az (zabs zr zi))
  32.       (setf iaz (f2cl-lib:int az))
  33.       (setf ifnu (f2cl-lib:int fnu))
  34.       (setf inu (f2cl-lib:int-sub (f2cl-lib:int-add ifnu n) 1))
  35.       (setf at (+ iaz 1.0))
  36.       (setf raz (/ 1.0 az))
  37.       (setf str (* zr raz))
  38.       (setf sti (* (- zi) raz))
  39.       (setf ckr (* str at raz))
  40.       (setf cki (* sti at raz))
  41.       (setf rzr (* (+ str str) raz))
  42.       (setf rzi (* (+ sti sti) raz))
  43.       (setf p1r zeror)
  44.       (setf p1i zeroi)
  45.       (setf p2r coner)
  46.       (setf p2i conei)
  47.       (setf ack (* (+ at 1.0) raz))
  48.       (setf rho (+ ack (f2cl-lib:fsqrt (- (* ack ack) 1.0))))
  49.       (setf rho2 (* rho rho))
  50.       (setf tst (/ (+ rho2 rho2) (* (- rho2 1.0) (- rho 1.0))))
  51.       (setf tst (/ tst tol))
  52.       (setf ak at)
  53.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  54.                     ((> i 80) nil)
  55.         (tagbody
  56.           (setf ptr p2r)
  57.           (setf pti p2i)
  58.           (setf p2r (- p1r (- (* ckr ptr) (* cki pti))))
  59.           (setf p2i (- p1i (+ (* cki ptr) (* ckr pti))))
  60.           (setf p1r ptr)
  61.           (setf p1i pti)
  62.           (setf ckr (+ ckr rzr))
  63.           (setf cki (+ cki rzi))
  64.           (setf ap (zabs p2r p2i))
  65.           (if (> ap (* tst ak ak)) (go label20))
  66.           (setf ak (+ ak 1.0))
  67.          label10))
  68.       (go label110)
  69.      label20
  70.       (setf i (f2cl-lib:int-add i 1))
  71.       (setf k 0)
  72.       (if (< inu iaz) (go label40))
  73.       (setf p1r zeror)
  74.       (setf p1i zeroi)
  75.       (setf p2r coner)
  76.       (setf p2i conei)
  77.       (setf at (+ inu 1.0))
  78.       (setf str (* zr raz))
  79.       (setf sti (* (- zi) raz))
  80.       (setf ckr (* str at raz))
  81.       (setf cki (* sti at raz))
  82.       (setf ack (* at raz))
  83.       (setf tst (f2cl-lib:fsqrt (/ ack tol)))
  84.       (setf itime 1)
  85.       (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
  86.                     ((> k 80) nil)
  87.         (tagbody
  88.           (setf ptr p2r)
  89.           (setf pti p2i)
  90.           (setf p2r (- p1r (- (* ckr ptr) (* cki pti))))
  91.           (setf p2i (- p1i (+ (* ckr pti) (* cki ptr))))
  92.           (setf p1r ptr)
  93.           (setf p1i pti)
  94.           (setf ckr (+ ckr rzr))
  95.           (setf cki (+ cki rzi))
  96.           (setf ap (zabs p2r p2i))
  97.           (if (< ap tst) (go label30))
  98.           (if (= itime 2) (go label40))
  99.           (setf ack (zabs ckr cki))
  100.           (setf flam (+ ack (f2cl-lib:fsqrt (- (* ack ack) 1.0))))
  101.           (setf fkap (/ ap (zabs p1r p1i)))
  102.           (setf rho (min flam fkap))
  103.           (setf tst (* tst (f2cl-lib:fsqrt (/ rho (- (* rho rho) 1.0)))))
  104.           (setf itime 2)
  105.          label30))
  106.       (go label110)
  107.      label40
  108.       (setf k (f2cl-lib:int-add k 1))
  109.       (setf kk
  110.               (max (the f2cl-lib:integer4 (f2cl-lib:int-add i iaz))
  111.                    (the f2cl-lib:integer4 (f2cl-lib:int-add k inu))))
  112.       (setf fkk (coerce (the f2cl-lib:integer4 kk) 'double-float))
  113.       (setf p1r zeror)
  114.       (setf p1i zeroi)
  115.       (setf p2r scle)
  116.       (setf p2i zeroi)
  117.       (setf fnf (- fnu ifnu))
  118.       (setf tfnf (+ fnf fnf))
  119.       (setf bk
  120.               (-
  121.                (multiple-value-bind
  122.                    (ret-val var-0 var-1)
  123.                    (dgamln (+ fkk tfnf 1.0) idum)
  124.                  (declare (ignore var-0))
  125.                  (setf idum var-1)
  126.                  ret-val)
  127.                (multiple-value-bind
  128.                    (ret-val var-0 var-1)
  129.                    (dgamln (+ fkk 1.0) idum)
  130.                  (declare (ignore var-0))
  131.                  (setf idum var-1)
  132.                  ret-val)
  133.                (multiple-value-bind
  134.                    (ret-val var-0 var-1)
  135.                    (dgamln (+ tfnf 1.0) idum)
  136.                  (declare (ignore var-0))
  137.                  (setf idum var-1)
  138.                  ret-val)))
  139.       (setf bk (exp bk))
  140.       (setf sumr zeror)
  141.       (setf sumi zeroi)
  142.       (setf km (f2cl-lib:int-sub kk inu))
  143.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  144.                     ((> i km) nil)
  145.         (tagbody
  146.           (setf ptr p2r)
  147.           (setf pti p2i)
  148.           (setf p2r (+ p1r (* (+ fkk fnf) (- (* rzr ptr) (* rzi pti)))))
  149.           (setf p2i (+ p1i (* (+ fkk fnf) (+ (* rzi ptr) (* rzr pti)))))
  150.           (setf p1r ptr)
  151.           (setf p1i pti)
  152.           (setf ak (+ 1.0 (/ (- tfnf) (+ fkk tfnf))))
  153.           (setf ack (* bk ak))
  154.           (setf sumr (+ sumr (* (+ ack bk) p1r)))
  155.           (setf sumi (+ sumi (* (+ ack bk) p1i)))
  156.           (setf bk ack)
  157.           (setf fkk (- fkk 1.0))
  158.          label50))
  159.       (f2cl-lib:fset (f2cl-lib:fref yr (n) ((1 n))) p2r)
  160.       (f2cl-lib:fset (f2cl-lib:fref yi (n) ((1 n))) p2i)
  161.       (if (= n 1) (go label70))
  162.       (f2cl-lib:fdo (i 2 (f2cl-lib:int-add i 1))
  163.                     ((> i n) nil)
  164.         (tagbody
  165.           (setf ptr p2r)
  166.           (setf pti p2i)
  167.           (setf p2r (+ p1r (* (+ fkk fnf) (- (* rzr ptr) (* rzi pti)))))
  168.           (setf p2i (+ p1i (* (+ fkk fnf) (+ (* rzi ptr) (* rzr pti)))))
  169.           (setf p1r ptr)
  170.           (setf p1i pti)
  171.           (setf ak (+ 1.0 (/ (- tfnf) (+ fkk tfnf))))
  172.           (setf ack (* bk ak))
  173.           (setf sumr (+ sumr (* (+ ack bk) p1r)))
  174.           (setf sumi (+ sumi (* (+ ack bk) p1i)))
  175.           (setf bk ack)
  176.           (setf fkk (- fkk 1.0))
  177.           (setf m (f2cl-lib:int-add (f2cl-lib:int-sub n i) 1))
  178.           (f2cl-lib:fset (f2cl-lib:fref yr (m) ((1 n))) p2r)
  179.           (f2cl-lib:fset (f2cl-lib:fref yi (m) ((1 n))) p2i)
  180.          label60))
  181.      label70
  182.       (if (<= ifnu 0) (go label90))
  183.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  184.                     ((> i ifnu) nil)
  185.         (tagbody
  186.           (setf ptr p2r)
  187.           (setf pti p2i)
  188.           (setf p2r (+ p1r (* (+ fkk fnf) (- (* rzr ptr) (* rzi pti)))))
  189.           (setf p2i (+ p1i (* (+ fkk fnf) (+ (* rzr pti) (* rzi ptr)))))
  190.           (setf p1r ptr)
  191.           (setf p1i pti)
  192.           (setf ak (+ 1.0 (/ (- tfnf) (+ fkk tfnf))))
  193.           (setf ack (* bk ak))
  194.           (setf sumr (+ sumr (* (+ ack bk) p1r)))
  195.           (setf sumi (+ sumi (* (+ ack bk) p1i)))
  196.           (setf bk ack)
  197.           (setf fkk (- fkk 1.0))
  198.          label80))
  199.      label90
  200.       (setf ptr zr)
  201.       (setf pti zi)
  202.       (if (= kode 2) (setf ptr zeror))
  203.       (multiple-value-bind
  204.           (var-0 var-1 var-2 var-3 var-4)
  205.           (zlog rzr rzi str sti idum)
  206.         (declare (ignore var-0 var-1))
  207.         (setf str var-2)
  208.         (setf sti var-3)
  209.         (setf idum var-4))
  210.       (setf p1r (+ (* (- fnf) str) ptr))
  211.       (setf p1i (+ (* (- fnf) sti) pti))
  212.       (setf ap
  213.               (multiple-value-bind
  214.                   (ret-val var-0 var-1)
  215.                   (dgamln (+ 1.0 fnf) idum)
  216.                 (declare (ignore var-0))
  217.                 (setf idum var-1)
  218.                 ret-val))
  219.       (setf ptr (- p1r ap))
  220.       (setf pti p1i)
  221.       (setf p2r (+ p2r sumr))
  222.       (setf p2i (+ p2i sumi))
  223.       (setf ap (zabs p2r p2i))
  224.       (setf p1r (/ 1.0 ap))
  225.       (multiple-value-bind
  226.           (var-0 var-1 var-2 var-3)
  227.           (zexp ptr pti str sti)
  228.         (declare (ignore var-0 var-1))
  229.         (setf str var-2)
  230.         (setf sti var-3))
  231.       (setf ckr (* str p1r))
  232.       (setf cki (* sti p1r))
  233.       (setf ptr (* p2r p1r))
  234.       (setf pti (* (- p2i) p1r))
  235.       (multiple-value-bind
  236.           (var-0 var-1 var-2 var-3 var-4 var-5)
  237.           (zmlt ckr cki ptr pti cnormr cnormi)
  238.         (declare (ignore var-0 var-1 var-2 var-3))
  239.         (setf cnormr var-4)
  240.         (setf cnormi var-5))
  241.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  242.                     ((> i n) nil)
  243.         (tagbody
  244.           (setf str
  245.                   (- (* (f2cl-lib:fref yr (i) ((1 n))) cnormr)
  246.                      (* (f2cl-lib:fref yi (i) ((1 n))) cnormi)))
  247.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n)))
  248.                          (+ (* (f2cl-lib:fref yr (i) ((1 n))) cnormi)
  249.                             (* (f2cl-lib:fref yi (i) ((1 n))) cnormr)))
  250.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n))) str)
  251.          label100))
  252.       (go end_label)
  253.      label110
  254.       (setf nz -2)
  255.       (go end_label)
  256.      end_label
  257.       (return (values nil nil nil nil nil nil nil nz nil)))))
  258.  
  259.